
#pragma once

#include "tbb/tick_count.h"
#include "tbb/task_scheduler_init.h"
#include "tbb/parallel_for.h"
#include "tbb/parallel_reduce.h"
#include "tbb/blocked_range.h" 

typedef Point point_t;
typedef tbb::blocked_range<size_t> range_t; 

#include "tbb/concurrent_vector.h"

typedef tbb::concurrent_vector<point_t> pointVec_t;

void appendVector(const point_t* src, size_t srcSize, pointVec_t& dest) {
	std::copy(src, src + srcSize, dest.begin() + dest.grow_by(srcSize));
}

void appendVector(const pointVec_t& src, pointVec_t& dest) {
	std::copy(src.begin(), src.end(), dest.begin() + dest.grow_by(src.size()));
} 

namespace cfg {
    // convex hull problem parameter defaults
    const long    NP = 10000000;			// problem size
    const int     SNT = 1;			// minimal number of threads
    const int     ENT = 8;			// maximal number of threads

    // convex hull problem user set parameters
    long   MAXPOINTS         = NP;
    int    NUM_THREADS_START = SNT;
    int    NUM_THREADS_END   = ENT;

    // convex hull grain sizes for 3 subproblems. Be sure 16*GS < 512Kb
    const size_t GENERATE_GS = 2500;
    const size_t FINDEXT_GS  = 2500;
    const size_t DIVIDE_GS   = 2500;
}; 

class FindXExtremum {
public:
    typedef enum {
        minX, maxX
    } extremumType;

    static const size_t  grainSize = cfg::FINDEXT_GS;

    FindXExtremum(const pointVec_t& points_, extremumType exType_)
        : points(points_), exType(exType_), extrXPoint(points[0]) {}

    FindXExtremum(const FindXExtremum& fxex, tbb::split)
        : points(fxex.points), exType(fxex.exType), extrXPoint(fxex.extrXPoint) {}

    void operator()(const range_t& range) {
        const size_t i_end = range.end();
        if(!range.empty()) {
            for(size_t i = range.begin(); i != i_end; ++i) {
                if(closerToExtremum(points[i])) {
                    extrXPoint = points[i];
                }
            }
        }
    }

    void join(const FindXExtremum &rhs) {
        if(closerToExtremum(rhs.extrXPoint)) {
            extrXPoint = rhs.extrXPoint;
        }
    }

    point_t extremeXPoint() {
        return extrXPoint;
    }

private:
    const pointVec_t    &points;
    const extremumType   exType;
    point_t              extrXPoint;
    bool closerToExtremum(const point_t &p) const {
        switch(exType) {
        case minX:
            return p.x<extrXPoint.x; break;
        case maxX:
            return p.x>extrXPoint.x; break;
        }
        return false; // avoid warning
    }
};

template <FindXExtremum::extremumType type>
point_t extremum(const pointVec_t &P) {
    FindXExtremum fxBody(P, type);
    tbb::parallel_reduce(range_t(0, P.size(), FindXExtremum::grainSize), fxBody);
    return fxBody.extremeXPoint();
}

class SplitByCP {
    const pointVec_t    &initialSet;
    pointVec_t          &reducedSet;
    point_t              p1, p2;
    point_t              farPoint;
    double               howFar;
public:
    static const size_t grainSize = cfg::DIVIDE_GS;

    SplitByCP( point_t _p1, point_t _p2,
        const pointVec_t &_initialSet, pointVec_t &_reducedSet)
        : p1(_p1), p2(_p2),
        initialSet(_initialSet), reducedSet(_reducedSet),
        howFar(0), farPoint(p1) {
    }

    SplitByCP( SplitByCP& sbcp, tbb::split )
        : p1(sbcp.p1), p2(sbcp.p2),
        initialSet(sbcp.initialSet), reducedSet(sbcp.reducedSet),
        howFar(0), farPoint(p1) {}

    void operator()( const range_t& range ) {
        const size_t i_end = range.end();
        double cp;
        for(size_t i=range.begin(); i!=i_end; ++i) {
            if( (initialSet[i] != p1) && (initialSet[i] != p2) ) {
                cp = cross_product(p1, p2, initialSet[i]);
                if(cp>0) { 
                    reducedSet.push_back(initialSet[i]); 
                    if(cp>howFar) {
                        farPoint = initialSet[i];
                        howFar   = cp;
                    }
                }
            }
        }
    }

    void join(const SplitByCP& rhs) {
        if(rhs.howFar>howFar) {
            howFar   = rhs.howFar;
            farPoint = rhs.farPoint;
        }
    }

    point_t farthestPoint() const {
        return farPoint;
    }
}; 
template <typename BodyType>
point_t divide(const pointVec_t &P, pointVec_t &P_reduced,
              const point_t &p1, const point_t &p2) {
    BodyType body(p1, p2, P, P_reduced);
    tbb::parallel_reduce(range_t(0, P.size(), BodyType::grainSize), body);
    return body.farthestPoint();
}

void divide_and_conquer(const pointVec_t &P, pointVec_t &H,
                        point_t p1, point_t p2) {
    if (P.size()<2) {
        H.push_back(p1); 
        appendVector(P, H); 
    }
    else {
        pointVec_t P_reduced;
        pointVec_t H1, H2;
        point_t p_far;

        p_far = divide<SplitByCP>(P, P_reduced, p1, p2);

        divide_and_conquer(P_reduced, H1, p1, p_far);
        divide_and_conquer(P_reduced, H2, p_far, p2);

        appendVector(H1, H);
        appendVector(H2, H);
    }
}

double tbbConvexHull(const Points &pts, Points &res)
{
	pointVec_t hull;
	pointVec_t points(pts.size());
	for (size_t i = 0; i < pts.size(); i++)
		points[i] = pts[i];
	tbb::tick_count t_start = tbb::tick_count::now();
    hull.clear();

    point_t p_maxx = extremum<FindXExtremum::maxX>(points);
    point_t p_minx = extremum<FindXExtremum::minX>(points);

    pointVec_t H;

    divide_and_conquer(points, hull, p_maxx, p_minx);
    divide_and_conquer(points, H, p_minx, p_maxx);
    appendVector(H, hull);
	tbb::tick_count t_end = tbb::tick_count::now();
	res.resize(hull.size());
	for (size_t i = 0; i < hull.size(); i++)
		res[i] = hull[i];
	return (t_end - t_start).seconds();
} 
